#!/usr/bin/env bash

##########################################################################################################################
## SCRIPT TO CALCULATE VOXEL-MIRRORED HOMOTOPIC CONNECTIVITY 
##
## Written by Xi-Nian Zuo, Maarten Mennes & Michael Milham 
## Edited by the Chief Underpants Gnome 12th April 2011
## for more information see www.nitrc.org/projects/fcon_1000
##
##########################################################################################################################


# USAGE

function usage {
    echo
    tput setf 3
    echo "Voxel-Mirrored Homotopic Connectivity: Single-Subject Processing"    
    echo    
    tput setf 2
    echo "Usage:"
    echo "  $0 <data dir> <subject list> <resting dirname>  <rest scan name> <registration type? LINEAR/NONLINEAR>"
    echo
    echo  "  EXAMPLE: ./singlesubjectVMHC.sh  /data/studyX  /data/studyX/scripts/subs.txt rest_1 rest YES NONLINEAR"     
    echo    
    tput setf 3   
    echo "Inputs:"         
    echo "  <data dir>: directory containing your preprocessed subject data (e.g., /data/studyX)"
    echo "  <subject list>: path to .txt file containing subject IDs (e.g., /data/studyX/scripts/subs.txt)"
    echo "  <resting dirname>: name of subject-level directory containing resting scan (e.g., rest_1)" 
    echo "  <rest scan name>: name of resting scan (e.g., rest)"    
    echo "  <registration type?>: Run linear or nonlinear registration to symmetric template? Enter LINEAR or NONLINEAR "
    echo 
    tput setf 7
}

if [[ "$#" = 0 ]]; then
    usage
    exit 2
fi


## analysisdirectory
dir=$1
## full/path/to/subject_list.txt containing subjects you want to run, you can also specify one subject
subject_list=$2
## func
func=$3
## name of resting-state scan (no extenstion)
rest=$4
##Registration type?
reg_type=$5

## standard brain as reference for registration
symm_standard_brain=${FSLDIR}/data/standard/MNI152_T1_2mm_brain_symmetric.nii.gz
symm_standard=${FSLDIR}/data/standard/MNI152_T1_2mm_symmetric.nii.gz


##########################################################################################################################
##---START OF SCRIPT----------------------------------------------------------------------------------------------------##
##########################################################################################################################

echo ---------------------------------------
echo !!!! CALCULATING VMHC !!!!
echo ---------------------------------------

cwd=$( pwd )


## Get subjects to run
if [ -f ${subject_list} ]
	then
	subjects=$( cat ${subject_list} )
	else
	echo "Subject list ${subject_list} does not exist"
fi


## A. SUBJECT LOOP
for subject in $subjects
	do

	image2use=${rest}_res_noPWf

	## directory setup

	func_dir=${dir}/${subject}/${func}
	reg_dir=${dir}/${subject}/${func}/2_reg
	anat_dir=${dir}/${subject}/anat
	anat_reg_dir=${anat_dir}/reg
	vmhc_dir=${dir}/${subject}/${func}/vmhc


	echo --------------------------
	echo running subject ${subject}
	echo --------------------------

	## check if func_dir is present
	if [ ! -d ${func_dir} ]
		then
		echo SUBJECT ${subject} has no ${func}
		continue
	fi

	## check if vmhc is already present, if not run subject
	if [ -f ${vmhc_dir}/VMHC.nii.gz -a -f ${vmhc_dir}/VMHC_Z.nii.gz -a -f ${vmhc_dir}/VMHC_Z_stat.nii.gz ]
		then
		echo SUBJECT ${subject} already has VMHC images
		echo SUBJECT ${subject} SKIPPED
		continue
	fi

	## make vmhc directory
	mkdir ${vmhc_dir}; cd ${vmhc_dir}

	## Run the registration to the symmetrical template
	if [ ! -f ${vmhc_dir}/${image2use}2symmstandard.nii.gz ]
		then

		echo using image registered to symmetrical standard
		echo starting registration
		
		if [ ${reg_type} = 'NONLINEAR' ]
			then
			
			if [ ! -f ${anat_reg_dir}/highres2symmstandard_warp.nii.gz ]
				then	
		
				## Linear registration of T1 --> symmetric standard
				echo Flirting T1 to symmetrical standard brain
				flirt -ref ${symm_standard_brain} -in ${anat_dir}/brain.nii.gz -out ${anat_reg_dir}/highres2symmstandard.nii.gz -omat ${anat_reg_dir}/highres2symmstandard.mat -cost corratio -searchcost corratio -dof 12 -interp trilinear

				## Perform nonlinear registration (higres to standard) to symmetric standard brain
				echo Fnirting T1 to symmetrical standard brain
				fnirt --in=${anat_dir}/head.nii.gz --aff=${anat_reg_dir}/highres2symmstandard.mat --cout=${anat_reg_dir}/highres2symmstandard_warp.nii.gz --iout=${anat_reg_dir}/fnirt_highres2symmstandard --jout=${anat_reg_dir}/highres2symmstandard_jac --config=T1_2_MNI152_2mm_symmetric --ref=${symm_standard} --refmask=${FSLDIR}/data/standard/MNI152_T1_2mm_brain_mask_symmetric_dil.nii.gz --warpres=10,10,10
			fi

			## Apply nonlinear registration (func to standard)
			echo Applying warp to ${image2use}
			applywarp --ref=${symm_standard} --in=${func_dir}/${image2use}.nii.gz --out=${vmhc_dir}/${image2use}2symmstandard.nii.gz --warp=${anat_reg_dir}/highres2symmstandard_warp.nii.gz --premat=${reg_dir}/example_func2highres.mat
	
		else
		
			if [ ! -f ${anat_reg_dir}/highres2symmstandard.nii.gz ]
				then	
				## Linear registration of T1 --> symmetric standard
				echo Flirting T1 to symmetrical standard brain
				flirt -ref ${symm_standard_brain} -in ${anat_dir}/brain.nii.gz -out ${anat_reg_dir}/highres2symmstandard.nii.gz -omat ${anat_reg_dir}/highres2symmstandard.mat -cost corratio -searchcost corratio -dof 12 -interp trilinear
			
				convert_xfm -omat ${reg_dir}/example_func2symmstandard.mat -concat ${anat_reg_dir}/highres2symmstandard.mat ${reg_dir}/example_func2highres.mat

				echo Applying linear registration to ${image2use}
				flirt -ref ${symm_standard_brain} -in ${func_dir}/${image2use}.nii.gz -out ${vmhc_dir}/${image2use}2symmstandard.nii.gz -applyxfm -init ${reg_dir}/example_func2symmstandard.mat
			fi
	
		fi
	fi


	## copy and L/R swap file
	echo swapping file
	fslswapdim ${vmhc_dir}/${image2use}2symmstandard.nii.gz -x y z tmp_LRflipped.nii.gz

	## caculate vmhc
	echo Calculating pearson correlation between res2standard and flipped res2standard
	3dTcorrelate -pearson -polort -1 -prefix VMHC.nii.gz ${vmhc_dir}/${image2use}2symmstandard.nii.gz tmp_LRflipped.nii.gz	
	
	## Fisher Z transform map
	3dcalc -a VMHC.nii.gz -expr 'log((1+a)/(1-a))/2' -prefix VMHC_Z.nii.gz
	## Z statistic map
	nvols=$( fslnvols tmp_LRflipped.nii.gz )
	echo ${nvols} volumes... 
	3dcalc -a VMHC_Z.nii.gz -expr 'a*sqrt('${nvols}'-3)' -prefix VMHC_Z_stat.nii.gz

	## Clean up
	echo cleaning up
	rm -rfv tmp*

done


